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Abstract 

We present and compare third- as well as fifth-order accurate finite difference schemes for the numerical solution of 
the compressible ideal MHD equations in multiple spatial dimensions. The selected methods lean on four different 
reconstruction techniques based on recently improved versions of the weighted essentially non-oscillatory (WENO) 
schemes, monotonicity preserving (MP) schemes as well as slope-limited polynomial reconstruction. The proposed 
numerical methods are highly accurate in smooth regions of the flow, avoid loss of accuracy in proximity of smooth 
extrema and provide sharp non-oscillatory transitions at discontinuities. 

We suggest a numerical formulation based on a cell-centered approach where all of the primary flow variables 
are discretized at the zone center. The divergence-free condition is enforced by augmenting the MHD equations 
with a generalized Lagrange multiplier yielding a mixed hyperbolic/parabolic correction, as in Dedner et al. (J. 
Comput. Phys. 175 (2002) 645-673). The resulting family of schemes is robust, cost-effective and straightforward to 
implement. Compared to previous existing approaches, it completely avoids the CPU intensive workload associated 
with an elliptic divergence cleaning step and the additional complexities required by staggered mesh algorithms. 

Extensive numerical testing demonstrate the robustness and reliability of the proposed framework for computa- 
tions involving both smooth and discontinuous features. 

Keywords: Magnetohydrodynamics, Compressible Flow, Higher-order methods, WENO schemes, Monotonicity 
Preserving, Cell-centered methods 



1. Introduction 

The development of high-order schemes has been receiving an increasing amount of attention from practitioners 
in the fields of fluid dynamics and, only more recently, magnetohydrodynamics (MHD). This interest is driven by a 
variety of reasons, such as the possibility of obtaining highly accurate solutions with reduced computational effort as 
well as the need to narrow the gap between the smallest resolved features and the dissipative scales. Although several 
successful strategies have been developed in the context of the Euler equations of gasdynamics, only few of them have 
been extended to MHD. In the present context, we focus our attention on high-order finite difference schemes for the 
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solution of the compressible MHD equations in multiple spatial dimensions, 

dp 



dt 



d( P y) 
dt 



+ V 



+ V • (pv) 



B 



pyy' -BB' + I \p + — 



8E V, 

<9f 



dB 
B 2 



V x (v x B) 



£ + /> + y|v-(v-B)B 



0. 
0. 
0. 



(1) 



where p, v, B, E and p are the fluid density, velocity vector, magnetic induction, energy and gas pressure, respectively. 
The system of equations ([D is complemented by the divergence-free constraint of the magnetic field, 



V-B = 0, 



(2) 



and by an equation of state relating energy and pressures. For the present work we assume an ideal gas law 

£ = f^T + 5( pv2 + B2 )' (3) 

where F is the ratio of specific heats. 

Traditional second-order schemes have been largely employed for the solution of Eq. 0} u sing either finite vol- 
ume (FV, e.g., BHmmUGllHliSS^r finite difference ( FD ' e -S-> ESElBlIiGlEi) methods. 
At the second-order level, the two approaches are essentially equivalent and popular schemes have been built on 
Godunov-type discretizations based on the Total Variation Diminishing (TVD, 12511 ) property making use of slope- 
limited reconstructions. In spite of the excellent results produced in proximity of discontinuous waves where sharp 
non-oscillatory transitions can be obtained, TVD schemes still suffer from excessive unwanted numerical dissipation 
in regions of smooth flow. This deficiency owes to the inherent behavior of TVD methods that reduces the order of 
accuracy to first-order near local extrema (clipping) and smear linearly degenerate fields (such as contact waves) much 
more than shocks. Furthermore, discretization errors are mainly responsible for the loss of accuracy. 

Efforts to relax the TVD condition and overcome these limitations have been spent over the last decades towards 
the development of highly accurate schemes that retain the robustness common to second-order Godunov-type meth- 
ods. The original piecewise parabolic method (PPM) method by lfl4ll . for example, provides fourth-order accurate 
interface values in smooth regions (in ID) and has been extended to MHD by ifl^flill and, more recently by J27, 281. 
PPM, however, still degenerates to first-order at smooth extrema and attempts to solve the problem have been recently 
presented in 111511 and 14511 . 

Based on a different approach, weighted essentially non-oscillatory (WENO, l47ll ) schemes have improved on their 
ENO predecessor (originally proposed by Harten et al. l26ll ) and are now considered a powerful and effective tool 
for solving hyperbolic partial differential equations. WENO methods provide highly accurate solutions in regions of 
smooth flow and non-oscillatory transitions in presence of discontinuous waves by combining different interpolation 
stencils of order r into a weighted average of order 2r - 1 . The nonlinear weights are adjusted by the local smoothness 
of the solution so that essentially zero weights are given to non smooth stencils while optimal weights are prescribed 
in smooth regions. WENO scheme have been formulated in the context of MHD using both FD Hll and FV 



formulations, 

s shush. 

Third- and fifth-order WENO schemes have been recently improved in terms of 
reduced dissipation, better resolution properties and faster convergence rates (see f53ll and |@1) and will be considered 
here. 

An alternative strategy is followed by the Monotonicity Preserving (MP) family of schemes by Suresh & Huynh 
1 48[] who proposed to carry the reconstruction step by first computing an accurate and stable interface value and then 
by imposing monotonicity- and accuracy-preserving constraints to limit the original value. MP schemes have been 
successfully merged with WENO methods by |4] and employed in the context of relativistic MHD by OS]. 

Finally, a reconstruction procedure that avoids the clipping phenomenon has been recently discussed by Cada & 
Torrilhon ifioll who devised a new class of nonlinear limiter functions based upon a non-polynomial reconstruction 
showing good shape-preserving properties. 
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It is important to point out that, for spatial accuracy higher than two, multidimensional FV schemes become 
notoriously more elaborate than their FD counterparts, since point values can no longer be interchanged with volume 
averages. As a result, FV schemes generally require fully multidimensional reconstructions and the solution of several 
Riemann problems at a zone face prov iding the necessary number of quadrature points required by the desired level 
of accuracy, see, for instance, Ill2ll49l Dll. However, FV algorithms do have the adventage that they are better suited 
to non-uniform grids and adaptive mesh hierarchies. High-order FV schemes have been recently ameliorated in the 
work of i2ll l7[l using either ADER-WENO schemes or least-squares polynomial reconstruction. 

Conversely, multidimensional FD schemes evolve the point values of the conserved quantities and considerably 
ease up the coding efforts by restricting the computations of flux derivatives to one dimensional stencils. In this 
perspective, we present a new class of FD numerical schemes adopting a point-wise, cell-centered formulation of 
all of the flow quantities, including magnetic fields. The proposed schemes have order of accuracy three and five 
and their performance is compared through extensive testing on two and three-dimensional problems. Selected third- 
order accurate schemes are i) an improved version of the classical third-order WENO scheme of I29I1 based on new 
weight functions designed to improve accuracy near critical points lf53ll and ii) the recentlyproposed non-polynomial 
reconstruction of 111 Oil . Selected fifth-order schemes include i) the WENO-Z scheme of 191] and ii) the monotonicity 
preserving scheme of 114 811 based on a fifth-order accurate interface value (MP5 henceforth). 

The solenoidal constraint of the mag netic field is controlled by extending the hyperbolic/parabolic divergence 
cleaning technique of Dedner et al. 12011 to FD schemes. This avoids the computational cost associated with an 
elliptic cleaning step as in J30I1 . and the scrupulous treatment of staggered fields demanded by constrained transport 



algorithms, e.g. |5|, |35l |27|, |7|] . Furthermore, Mignone & Tzeferacos 1 39] have shown through extensive testing, for a 
class of second-order accurate schemes, that the GLM approach is robust and can achieve accuracy comparable to the 
constrained transport. The resulting class of schemes is explicit and fully conservative in mass, momentum, magnetic 
induction and energy. Besides the ease of implementation and efficiency issues, the benefits offered by a method 
where all of the primary flow variables are placed at the same spatial position ease the task to add more complex 
physics. 

The comparison between the different methods of solution is conveniently handled using the PLUTO code for 
computational astrophysics 13711 . 

The paper is structured as follows. In §f2] we describe the GLM-MHD equations, while ^3] shows the finite dif- 
ference formulation and the selected reconstruction methods. In §|4]we test and compare the different scheme perfor- 
mance on problems involving the propagation of both continuous and discontinuous features. Conclusions are drawn 
in SSI 



2. The Constrained GLM-MHD Equations 

We look at a conservative discretization of the MHD equations ([]]) where all fluid variables retain a cell-centered 
collocation and enforce the divergence-free condition through the hyperbolic/parabolic divergence cleaning technique 
of Dedner's lEoll . In this approach Gauss's and Faraday's laws of magnetism are modified by the introduction of a new 
scalar field function or generalized Lagrangian multiplier (GLM henceforth) if/. The resulting system of GLM-MHD 
equations then reads 

<9U v d¥, c 

l=x,y,z 

with conservative state vector U and fluxes F/ defined by 
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where d = x,y,z labels the different components while Sdi is the delta Kronecker symbol. Equations © are hyperbolic 
and fully conservative with the only exception of the unphysical scalar field \p which satisfies a non-homogeneous 
equation with a source term. In the GLM approach, divergence errors are propagated to the domain boundaries at 
finite speed q, and damped at a rate given by cjjc^ (see 

The eigenvalues of the MHD flux Jacobians (9F//3U are all real and coincide with the ordinary MHD waves plus 
two additional modes ±Ch, for a total of 9 characteristic waves. Restricting our attention to the I = x direction, they 
are given by 



A 1,9 = +c h , A 2 '* = v x + c f , A X1 = v x + c a , A 4 ' 6 = v x + c s , A 5 = v x , (6) 



where 



are the fast magneto-sonic (c/ with the + sign), slow magneto-sonic (c s with the - sign) and Alfven velocities. The 
two additional modes +c;, are decoupled from the remaining ones and corresponds to linear waves carrying jumps in 
B x and if/. These waves are made to propagate at the maximum signal speed compatible with the time step, i.e., 

c h = max (\v x \ + Cfe, \v y \ + c/, v , \v z \ + c f ~) . (8) 

where Cf sX ,Cfy,Cf iZ are the fast magneto-sonic speeds in the three directions and the maximum is taken throughout the 
domain. 

Owing to the decoupling, one can treat the 2x2 linear system given by the longitudinal component of the field B/ 
and tp separately from the other ordinary 7-wave MHD equations. As we shall see, this greatly simplifies the solution 
process and allows to use the standard characteristic decomposition of the MHD equations. 

Following ll39ll . we divide the solution process into an homogeneous step, where the GLM-MHD © are solved 
with S = 0, and a source step, where integration is done analytically: 

r ,(Af) _ ,,,((>>-. c h \ „ , Ch 



iffW =^)^ v \- ap —SL- \ , with a p = Ah-£. (9) 



- p 



where Ah = min(Ajc, Ay, Az) is the minimum grid size. Extensive numerical testing has shown that divergence errors 
are minimized when the parameter a p lies in the range [0,1] depending on the particular problem, although in presence 
of smooth flows this choice seems to be less sensitive to the numerical value of a p . 



3. Finite Difference schemes 

We consider a conservative finite difference discretization of (|4]i where point-values rather than volume averages 
are evolved in time. A uniform Cartesian mesh is employed with cell sizes Ax x Ay x Az centered at (xi,yj,Zk), 
where k label the computational zones in the three directions. For clarity of exposition, we disregard the integer 
subscripts when redundant but always keep the half increment index notation when referring to a cell boundary, e.g., 

Integration in time resorts to a semi-discrete formulation where, given a high-order numerical approximation -C(U) 
to the derivatives appearing on the right hand side of Eq. one is faced with the solution of the following initial 
value problem 

^=X(U), (10) 
at 

with initial condition given by the point-wise values of U(jc,-, yj,Zk, t n ) = V" - k - We choose the popular third-order 
Runge-Kutta scheme 1 46[ 24 ] to advance the solution in time, for which one lias 



U* = U"+£(U") 

444 



U" = T U" + T IT + — £(IT) , 



U" +1 = + j\J"* + jAt"£(\]") . 



The choice of the time step At" is restricted by the Courant-Friedrichs-Levy (CFL) condition: 



Af"=C fl — , (12) 

Ch 

where C a is the CFL number. Since the time step is proportional to the mesh size, the overall accuracy of the scheme 
is restricted to third-order because of the time-stepping introduced in Eq. (fTTT i. 

Our task is now to provide a stable and accurate non-oscillatory numerical approximation to -C(U). To this purpose, 
we begin by focusing our attention to the x- direction and set, for ease of notations, F,- = V x (Uijjc)- We then let point 
values of the flux F, correspond to the volume averages of another function, say F, and define 

Fi = hP +if( ® d t = ~L\ RiXi+ l ) ~ 11iXi -l ) \ , where HW = J t( ^- (13) 

In this formalism, point values of the flux F, are identified as cell averages of F(x) and H(x) may be regarded as the 
primitive function of F. Straightforward differentiation of Eq. ( fT3l yields the conservative approximation 



dx 



= s(*»i-*h)- ^ 



Stated in this form, the problem consists of finding a high-order approximation to the interface values of F j+ 1 knowing 
the undivided differences of the primitive function H(x), a procedure entirely analogous to that used in the context of 
finite volume methods such as PPM fl4ll . Thus one can set 

F i+ , =K(F W ), (15) 

where 7?() is a highly accurate reconstruction scheme providing a stable interface flux value from point-wise values 
and the index [s] spans through the interpolation stencil. 

The procedure can be repeated in an entirely similar way also for the y and z flux contributions and allows to write 
the X. operator in dTOb as 

-ecu) = ~ («u 4 - f„ .) - ^ (f,„ - Kj-i) - ^ (f z , h - Kh) ■ ^ 

This yields the fully unsplit approach considered in this paper. Alternatively, one could use a directionally split 
formalism to obtain the solution through a sequence of one dimensional problems separately corresponding to each 
term in equation 1161 . 

In order to ensure robustness and to avoid the appearance of spurious oscillations, the reconstruction step is best 
carried with the help of local characteristic fields and by separately evaluating contributions coming from right- and 
left-going waves. To this end we first compute, using the simple arithmetic average U i+ i = (U,- + U, + i)/2, left and 
right eigenvectors I/ + , and R* + , of the Jacobian matrix <9F/<9U, for each characteristic field k = 1, . . . , 9. We then 

obtain a projection of the positive and negative part of the flux using a simple Rusanov Lax-Friedrichs flux splitting: 



\U. .(F w +oTJ w ) 



(17) 



where F[ S ] and U[ S ] are the point-wise values of the flux and conservative variables. For a typical one-point upwind- 
biased approximation of order (2r + 1), one has [s] = i - r, . . . , i + r while [s'] = 2i — [s] + 1 mirrors left-going 
characteristic fields with respect to the interface i + \- The coefficient a" represents the maximum absolute value of 
the K-th characteristic speed throughout the domain. 

The global Lax-Friedrichs flux splitting thus introduced is particularly diffusive and other forms of splitting are of 
course possible, e.g. 11291 1411. However, we have found that the level of extra numerical dissipation tend to become less 
important for higher-order scheme. 
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The interface flux is then written as a local expansion in the right-eigenvector space: 



where the coefficients 



K 



(18) 



(19) 



are the reconstructed interface values of the local characteristic fields and 7?() can be any one of the procedures 
described in < 



3.1. Modification for the Constrained GLM-MHD equations 

The procedure illustrated so far is valid for an arbitrary system of hyperbolic conservation laws, provided L* and 
R* satisfy 

<9F 

L* • — ■ R k = A , (20) 
ov 

i.e., they are left and right eigenvectors of the flux Jacobian, respectively. However, following 12011 . we wish to 
exploit the full 7x7 characteristic decomposition of the usual MHD equations rather than resorting to a full 9x9 
diagonalization procedure. To this purpose, we take advantage of the fact that the longitudinal component of the field 
B x and the Lagrange multiplier ip satisfy 
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and are thus decoupled from the remaining seven MHD equations. Eq. (l2"TT i defines a constant coefficient linear 
hyperbolic system with left and right eigenvectors given, respectively, by the rows and columns of 



(22) 



associated with the eigenvalues A 1 = -c/, and A 9 = +c/,. The 2x2 linear system d2"TT i can be preliminary solved to find 
the values of B x and if/ at a given interface. Indeed, by applying the projection (TTTl i to the linear system (l2"TT i using Eq. 
one obtains that the only non trivial characteristic fields are 
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Since the eigenvectors are constant in space, the local projection at i + | are completely unnecessary and the compu- 
tations in Eq. ( 1231 can be carried out very efficiently throughout the grid. Once d23l have been reconstructed using 
Eq. ( fT9l one defines 

B xH={%-vt~i) /Ch ' ^H=% + %' (24) 

and proceed by solving the ordinary 7x7 MHD equations using B K j+ i defined by (f24T > as a constant parameter. 

3.2. Third and Fifth-order Accurate Reconstructions 

We have shown in §[3] that flux derivatives may be written in conservative form by applying any one-dimensional 
finite volume reconstruction to the point values of the flux F, . Among the variety of different strategies we investigate 
both third- and fifth-order accurate interpolation schemes making use of three- and five-point stencil, respectively: 



• an improved version of the classical third-order WENO scheme of |29] based on new weight functions designed 
to improve accuracy near critical points (WENO+3, 33.2. U ; 

• the recently proposed Lim03 third-order reconstruction of ifioll . 33.2.21 
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the improved WEN05 scheme of H also known as WENO-Z ( KDl ; 



the monotonicity preserving scheme of 1 48] based on a fifth-order interface value (MP5, 33.2.41 ). 



Our choice is motivated by the sake of comparing well-known and recently presented state of the art algorithms that 
rely on heavy usage of conditional statements (Lim03 and MP5) or completely avoid them (WENO+3 and WENO-Z). 

The proposed algorithms are applied to the left (-) and right (+) propagating characteristic fields defined by Eq. 
( TTTb to provide an accurate interface value, formally represented by Eq. ( fT9l ). Thus, in our formulation, the total 
number of reconstruction is 16: two for the linear characteristic fields defined by Eq. d23l and 14 for the left- and 
right-going wave families defined by Eq. ( TTTb with k = 2, . . . , 8. 

In the following we will drop the i + \ index for the sake of exposition and shorten either one of dTvT > with fyy 
Undivided difference will be frequently used and denoted with 

A i+i =f M -fi. (25) 
Occasionally, we will also make use of the Minmod and Median functions defined, respectively as 

Minmod(a, b) = ~ ~ ~ mm (M> 1^1) > Median(a, b,c) — a + Minmod(Z? - a,c - a) . (26) 

3.2.1. Third-Order Improved WENO (WENO+3) 

In the classical third-order WENO scheme of J29ll . the interface value is reconstructed using the information 
available on a three-point local stencil (xi-\,Xi, Xj+i). More specifically, a third-order accurate value is provided by a 
linear convex combination of second-order fluxes: 

R Vis] ) = wo + oi\ . (27) 



The weights to/ for / = 0, 1 are defined by 



co, = , a, = - - x2 , with p a =A 2 lt fa = A 2 , , (28) 

where do = 2/3, d\ = 1/3 are optimal weights and the smoothness indicators fii give a measure of the regularity of the 
corresponding polynomial approximation. 



The scheme has been recently improved in the work by Yamaleev & Carpenter, 115 311 . where the introduction of an 
additional nonlinear artificial dissipation term was shown to make the scheme stable in the L2-energy norm for both 
continuous and discontinuous solutions. Yamaleev & Carpenter also derived new weight functions providing faster 
convergence and improved accuracy at critical points. The improved weights are still defined by Eq. ( 1281 with 07 
replaced by 



at 



1 



A i+ |-A ; _ 



(29) 



To avoid loss of accuracy at critical points, it was shown in [53] that e has to satisfy e = 0(Ax 2 ). 

Here adopt the conventional third-order scheme defined by Eq. (I2"7li-(l2"8ll but with 07 replaced by Eq. d29l and 



simply set e = Ax 2 . This improves the accuracy over the original 3 rd order scheme of B29Q in regions where the solution 
is smooth and provides essentially non-oscillatory solutions near strong discontinuities and unresolved features. The 
improved third-order WENO scheme just described will be referred to as WENO+3. 



3.2.2. Third-Order Limited reconstruction (Lim03 ) 

Recently, Cada and Torrilhon ifToll have proposed a new and efficient third-order limiter function in the context 
of finite volume schemes. Similarly to the 3 rd -order WENO scheme described in 33.2.11 the new limiter employs a 
local three-point stencil to achieve piecewise-parabolic reconstruction for smooth data and preserves the accuracy at 
local extrema, thus avoiding the well known clipping of classical second-order TVD limiters. Interface values are 
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reconstructed using a simple piecewise-linear max/min function acting as a logical switch depending on the left and 
right slope: 

A, 



ft (/m) = fi + -y- [pm +x{m - p 3 (0))] 



(30) 



where 9 = A,_i/A l+ i is the slope ratio, Pt>(6) — (2 + ff)/3 is the building block giving polynomial quadratic recon- 
struction and <h(ff) is the third-order limiter 



4,(6) 



max [0,min (P 3 (ff),20, 1.6)] 
0\ 



0, min 



if > , 
if 6» < . 



(31) 



The function x in Eq. (f30l > smoothly switches between limited and unlimited reconstructions based on a local 
indicator function r\ properly introduced to avoid loss of accuracy at smooth extrema with one vanishing lateral 
derivative: 

,\i A 2 , + A 2 

X]- 



X = max 



1 n- 1 
0,mn 1,- + 

2 2e 



(32) 



' (rA*) 2 ' 

where e = 10~ 12 . The function rj measures the curvature of non-monotone data inside a computational zone and the 
free-parameter < r < 1 is used to discriminate between smooth extrema and shallow gradients. Larger values of r 
noticeably improve the reconstruction properties at the cost of introducing more local variation, see 111 Oil . In the tests 
presented here we use r — 1 . 

3.2.3. Fifth-Order Improved WENO: WENO-Z 

Borges et al. |@] presented an improved version of the classical fifth-order weighted essentially non-oscillatory 
(WENO) FD scheme of i29ll . The new scheme, denoted with WENO-Z, has been shown to be less dissipative and 
provide better resolution at critical points at a very modest additional computational cost. We will employ such 
scheme here and, for the sake of completeness, report only the essential steps for its implementation (for a thorough 
discussion see the paper by 0). 

Following the general idea of WENO reconstruction, one considers the convex combination of different third-order 
accurate interface values built on the three possible sub-stencils of i - 2 < s < i + 2: 



2/i- 2 - 7/m + nfi , -f-x + 5ff + 2f i+l 



■ u>\ - 



+ C02 



2f t + 5f M -f M 



(33) 



The weights w; for / = 0, 1 , 2 are defined by 
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(34) 



where do = 1/10, di = 3/5,^2 = 3/10 are the optimal weights giving a fifth-order accurate approximation, e = 10~ 40 
is a small number preventing division by zero and the smoothness indicators /?/ give a measure of the regularity of the 
corresponding polynomial approximation: 
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(35) 



While maintaining the essentially non-oscillatory behavior, the new formulation makes use of higher-order informa- 
tion about the regularity of the solution thus providing enhanced order of convergence at critical points as well as 
reduced dissipation at discontinuities. 
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3.2.4. Fifth-order Monotonicity Preserving (MP5) 

The monotonicity preserving (MP) schemes of Suresh & Huynh 14811 achieve high-order interface reconstruction 
by first providing an accurate polynomial interpolation and then by limiting the resulting value so as to preserve 
monotonicity near discontinuities and accuracy in smooth regions. The MP algorithm is better sought on stencils with 
five or more points in order to distinguish between local extrema and a genuine 0(1) discontinuities. Here we employ 
the fifth-order accurate scheme based on the (unlimited) interface value given by 

, 2f i - 2 -l3f i - l +41f i + 21f i+l -3f i+2 

/»i go ' (36) 

based on the five point values /;_2, . . . , Together with d36b . we also define the monotonicity-preserving bound 

/ MP =fi + Minmod(A i+ i,arA i _i) , (37) 

resulting from the median between fj, fi+\ and the left-sided extrapolated upper limit / UL = f, ■ + aA t _i. The parameter 
a > 2 controls the maximum steepness of the left sided slope and preserves monotonicity during a single Runge-Kutta 
stage (Eq. fTTb provided the CFL number satisfies C a < 1/(1 + a). In practice, setting a = 4 still allows larger values 
of C a to be used. The interface value given by Eq. d36t is not altered when the data is sufficiently smooth or monotone 
that f t+ i lies inside the interval defined by [//,/ MP ]. Otherwise limiting takes place by bringing the original value 
back into a new interval /[/ mm ,/ inax ] specifically designed to preserve accuracy near smooth extrema and provide 
monotone profile close to discontinuous data. The final reconstruction can be written as 

[ f i+ i if (f i+ i-Mfi + i-f MP )<0, 

K(fls ] ) = \ . . \ (38) 

[ Median (f mm ,f i+ 1, / max ) otherwise, 

where 

r m = max [min (f h f M , / MD ) , min (/, , / UL , / LC )1 , 

(39) 

/ max = min [max (f h f M , f MD ) , max (/,, / UL , / LC )] . 

The bounds given by Eq. (l39b provide accuracy -preserving constraints by allowing the original interface value f t+ i to 
lie in a somewhat larger interval than or /[/i,/ UL ]. This is accomplished by considering the intersection of 

the two extended intervals I[fj, f, + \, / MD ] and /[/i, / UL , / LC ] that leave enough room to accommodate smooth extrema 
based on a measure of the local curvature defined by 

_ m inmod (4d ; - d i+ 1 , 4d i+ i-du d, , d i+ 1 ) , (40) 
where dj = A i+ 1 - A,_ i . Using Eq. d40l i. one defines the median / MD and the large curvature / LC values as 

/M p = fi±fM_ _ 1 M4 /L C =f 1a , i^M4 (41) 

respectively. The curvature measure provided by (1401 is somewhat heuristic and chosen to reduce the amount of 
room for local extrema to develop. The reconstruction illustrated preserves monotonicity and does not degenerate to 
first-order in proximity of smooth extrema. 



4. Numerical Tests 

In this section we present a series of test problems aimed at the verification of the FD methods previously de- 
scribed. The selected algorithms have been implemented in the PLUTO code for astrophysical gas-dynamics 0711 in 
order to ease inter-scheme comparisons through a flexible common computational framework. 
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Unless otherwise stated, the specific heat ratio will be set to T = 5/3 and the Courant number C a will be taken 
equal to 0.8, 0.4 or 0.3 for one, two and three dimensional computations, respectively. Errors for a generic flow 
quantity Q are computed using the L\ discrete norm defined by 



(42) 



i.j.k 



where the summation extends to all grid zones, N x , N y and N z are the number of grid points in the three directions and 
Q le{ is a reference solution. The divergence of magnetic field is quantified using Eq. ( l42l with V ■ B computed as 



V B 



B y,M~ B yJ- 



Ax Ay 
where the interface values are obtained through Eq. ( l24l i. 



Az 



(43) 



4.1. Propagation of Circularly polarized Alfven Waves 

We start by considering a planar, circularly polarized Alfven wave propagating along the x direction. As the wave 
propagates, density and pressure stay constant whereas transverse vector components trace circles without changing 
their magnitude. Denoting with <d and k the angular frequency and wavenumber, respectively, one has 



(44) 



where <p — kx - a>t, a>/k — vo x ± c a is the corresponding phase velocity (c a is the Alfven speed) and A is the wave 
amplitude. The plus or minus sign corresponds to right or left propagating waves, respectively. Here we consider a 
standing wave for which one has vq x = vo v = vq z — and further setp = 1, c a = 1. 

The one-dimensional solution given by fl4t is first rotated by an angle y around the y axis and subsequently by 
an angle a around the z axis, as in l39ll . The resulting transformation leaves scalar quantities invariant and produces 
vector rotations q — > R yff q, where q is either velocity or magnetic field and 



' v x - 




v * 




' B x ' 






Vy 




vov + A sin 




By 




+ yJpA sin <p 






v \'0z + A COS (f> , 




k B z ) 




K + y/pA COS (p , 



' cos a cos y - sin a - cos a sin y s 
sin a cosy cos a - sine sin y 
sin y cosy 



cos a cosy 

- sin a 
- cos a sin y 



sin a cosy siny ' 

cos a 
- sin a sin y cos y J 



(45) 



are the rotation matrix and its inverse. 

Note that the rotation can be equivalently specified by prescribing the orientation of the wave vector k = (k x , k y , k z ) 
in a three-dimensional Cartesian frame through the angles a and /? such that 



k y k- 
tan a = — , tan 8 — — , 

k x k x 



(46) 



such that tany = cosatan/3. With these choices, <p in d44l becomes (p = k -x - cot where to = ±|k|. 

Periodicity is guaranteed by setting, without loss of generality, k x = 2n and by choosing the computational domain 
x G [0, 1], y G [0, 1/ tano-] and z G [0, 1/ tan/?]. With these definitions the wave returns into the original position after 
one period 

T = . = = (47) 
yl + tan- a + tan 2 /? 

Different configurations can be specified in terms of the four parameters a,/3,A and po (background pressure). One 
and two dimensional propagation are recovered by setting a — /3 = and /3 = 0, respectively. 
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Table 1 : Accuracy analysis for the one dimensional (third and fourth columns) and three dimensional (fifth and sixth) Alfven wave propagation 



after one wave period. Errors are computed as (B A ) 2 + ei (B v ) 2 + £1 (B z ) 2 . The numerical scheme and the number of points in the x direction 
are given in the first and second columns. For 3-D propagation, the resolution in the y and z direction is set by N y = N. = N x /2. 







One Dimension 




Three Dimensions 




Method 


N x 


ei (|B|) 


Ou 


ex (|B|) 


Ou 


WENO+3 


16 


3.45E-03 




2.54E-02 






32 


4.39E-04 


2.97 


3.68E-03 


2.79 




64 


5.52E-05 


2.99 


4.47E-04 


3.04 




128 


6.91E-06 


3.00 


5.51E-05 


3.02 




256 


8.64E-07 


3.00 


6.85E-06 


3.01 


Lim03 


16 


3.36E-03 




2.82E-02 






32 


4.36E-04 


2.95 


3.76E-03 


2.91 




64 


5.53E-05 


2.98 


4.34E-04 


3.11 




128 


6.91E-06 


3.00 


5.46E-05 


2.99 




256 


8 65F-07 


3.00 


fi 84F-06 

VJ.OT^J— i \J\J 


3.00 


WENO-Z 


16 


7.50E-04 




4.10E-03 






32 


2.40E-05 


4.96 


1.32E-04 


4.96 




64 


7.55E-07 


4.99 


3.89E-06 


5.09 




128 


2.36E-08 


5.00 


1.20E-07 


5.02 




256 


7.37E-10 


5.00 


3.74E-09 


5.00 


MP5 


16 


7.38E-04 




3.41E-03 






32 


2.40E-05 


4.94 


1.19E-04 


4.84 




64 


7.55E-07 


4.99 


3.81E-06 


4.97 




128 


2.36E-08 


5.00 


1.20E-07 


4.99 




256 


7.37E-10 


5.00 


3.74E-09 


5.00 



4.1.1. One Dimensional Propagation 

As a first test, we consider one-dimensional propagating waves on the segment x e [0, 1] using N x = 2 q grid 
points with q = 4, . . . , 8. We set the background pressure to be po = 0.1 and the wave amplitude A = 0.1. The GLM 
correction is not necessary and has turned off for one dimensional propagation. 

In order to investigate the convergence of solution, the integration time step is adjusted to 

M N = At No {^j P (48) 

where Atm is the nominal time increment at the minimum resolution No, whereas r > 3 is the spatial accuracy of the 
scheme. Errors (in L\ norm) for the four selected schemes are plotted after one wave period T = 1 in the left panel 
of Fig I A. 1 1 and arranged, together with the corresponding order of convergence, in the third and fourth columns of 
TableQ] All schemes meet the expected order of accuracy (i.e. 3 for Lim03 and WENO+3, 5 for WENO-Z and MP5) 
with no significant differences. It is remarkable that, at the resolution of 64 zones, the fifth-order schemes achieve 
essentially the same accuracy as the third-order schemes that make use of four times (i.e. N x = 256) as many points. 

4.1.2. Three Dimensional Oblique Propagation 

A three dimensional configuration is obtained by rotating the one-dimensional setup described in {s!4.1.1l bv the 
angles a = fl = tan -1 2 so that tany = 2/ V5 in Eq. (1451 1. The background pressure is po = 0.1 and the wave has 
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amplitude A = 0.1. The size of the computational box turns out to be x e [0,1], y e [0, 1/2], z e [0, 1/2] and the 
number of grid points is set by N v = N z = N x /2, where N x changes as in 34.1.11 Integration lasts for one wave period, 
i.e., t = T = 1/3 and the time step is determined by the same condition given by Eq. (148b . Thus, apart from the 
different normalization, our setup is identical to that used in 12811 . 

Errors are plotted at different resolutions in the right panel of Fig IA. II and sorted in Table [TJ for all schemes. On 
average, errors are ~ 4 larger than their one-dimensional counterparts but the overall behavior meets the expected 
order of accuracy with MP5 and Lim03 performing slightly better than WENO-Z and WENO+3, respectively. As for 
the ID case, roughly 1 /4 of the resolution is required by a fifth-order scheme to match the accuracy of a third-order 
one. 



Following [28], we construct in Fig IA. 21 a scatter plot of the magnetic field component parallel to the y axis of 
the original one dimensional frame. This is achieved by plotting, for every point in the computational domain, the y 
component of Ry t 'B as a function of the normal (x) coordinate of R^x, where R yff is the rotation matrix introduced 
in (l45l l. The ability of the scheme to retain the planar symmetry during the computation is confirmed by the lack of 
scatter in the plots. The profiles at different resolutions verify the general trend established in Table[TJand deviations 
from the exact solution appear to be imperceptible for N x > 64 for the third-order schemes and already at > 32 for 
the fifth-order schemes. 

Overally, the results obtained with third- and fifth-order accurate schemes outperform traditional TVD schemes, 
such as the CT-PPM algorithm of [28] yielding at most second-order accurate solutions. The CPU costs associated 
with WENO+3, Lim03 , WENO-Z and MP5 show, for this test problem, a relative scaling 1 : 0.98 : 1.46 : 1.31, 
respectively. 



4.1.3. Numerical Dissipation and Long Term Decay in Two Dimensions 

As already stated, circularly polarized Alfven waves are an exact nonlinear solution of the MHD equations and 
measuring their decay provides a direct indication of the intrinsic numerical viscosity and resistivity possessed by the 
underlying algorithm, see Ji^ElS]- This study is relevant, for example, in the field of MHD turbulence modeling 
where one should carefully control the amount of directionally-biased dissipation introduced by waves propagating 
inclined to the mesh. The error introduced during an oblique propagation is usually minimized at 45° since contri- 
butions coming from different directions have comparable magnitude. On the contrary, waves propagating at smaller 
inclination angles make the problem more challenging. 

Our setup builds on |5[] although we adopt a slightly different, more severe, configuration. Using the notations 
introduced in 34.11 we set tana- = 6, tan/3 = 0, A — 0.2 and prescribe the background pressure to be po = 1. The 
corresponding ratio of the plasma pressure to the (unperturbed) magnetic pressure is then given by p/(2pcl) = 1/2, 
where c a — 1 is the wave propagation speed. The choice of the inclination angle determines the computational domain 
x e [0, 1], y e [0, 1/6] as well as the wave period T — 1/ V37 from Eq. ( 1471 ). The final integration time t — 16.5 is 
chosen by having the wave cross the domain » 100 times. This configuration results in a more arduous test than ^ 
where the wave period was 6 V4~7r longer and the integration was stopped after * 37 wave transits. 

Fig IA. 31 shows, at the resolution of 120 x 20 mesh points, the maximum values of the vertical z components of 
velocity (left panel) and magnetic field (right panel) as functions of time. By the end of the simulation, third-order 
schemes (dashed lines) show some degree of dissipation with the wave amplitude being reduced to ~ 20 per cent of 
its initial value. On the contrary, schemes of order five (solid lines in the figure) preserve the original shape more 
accurately and the amplitude retains ~ 94 per cent of its nominal value. These results are compared, for illustrative 
purposes, to a 2 nd order TVD scheme using the Monotonized Central difference limiter (dotted lines), showing that the 
initial peak values have scaled down to ~ 3 per cent, thus showing a considerably larger level of numerical dissipation. 

These results are in agreement with previous investigations Q 0] and strongly supports the idea that problems 
involving complex wave interactions may benefit from using higher-order schemes such as the ones presented here. 

4.2. Shock tube problems 

Shock tube problems are commonly used to test the ability of the scheme in describing both continuous and 
discontinuous flow features. In the following we consider two and three dimensional rotated configurations of standard 
one dimensional tubes. The default value for the parameter a p controlling monopole damping (see Eq. |9]l is 0.8. 
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Table 2: One dimensional L\ norm errors for density, normal component of magnetic field and |V ■ B| for the two and three-dimensional shock tube. 

Two Dimensions Three Dimensions 



Method 


ei (p) 


eiCBi) 


ex (V ■ B) 


ei(p) 


ei(Bi) 


ex (V ■ B) 


WENO+3 


4.11E-03 


8.53E-05 


7.19E-03 


1.82E-03 


4.41E-05 


7.12E-03 


Lim03 


3.61E-03 


8.74E-05 


1.08E-02 


1.63E-03 


4.07E-05 


9.41E-03 


WENO-Z 


2.72E-03 


7.90E-05 


1.48E-02 


1.29E-03 


5.41E-05 


1.59E-02 


MP5 


2.31E-03 


6.24E-05 


7.60E-03 


1.07E-03 


2.22E-05 


1.37E-02 



i.2.1. Two-Dimensional Shock tube 

Following [30[ 33], we consider a rotated version of the Brio-Wu test problem Ji 
'iven by 

V L = (1,0, 0,0.75, l,l) r for xi<0, 

V s = (0.125,0,0,0.75, -1,0.1/ for *i>0, 



with left and right states are 



(49) 



where V = (p, vi, V2, B\,B2,p) is the vector of primitive variables. The subscript "1" gives the direction perpendicular 
to the initial surface of discontinuity whereas "2" corresponds to the transverse direction. Here r = 2 is used and the 
evolution is interrupted at time t = 0.2, before the fast waves reach the borders. 

In order to address the ability to preserve the initial planar symmetry we rotate the initial condition by the angle 
a = 7r/4 in a two dimensional plane with x e [-1, 1] and y e [-0.01,0.01] using x A^/100 grid points, with 
N x = 600. Vectors follow the same transformation given by Eq. ( t45T > with B = y — 0. This is known to minimize errors 
of the longitudinal component of the magnetic field (see for example the discussions in J52 , 27 ]). Boundary conditions 
respect the translational invariance specified by the rotation: for each flow quantity we prescribe q(i, j) = q{i+5i, j+6 j) 
where (Si,5j) = (1,-1), with the plus (minus) sign for the leftmost and upper (rightmost and lower) boundary. 
Computations are stopped before the fast rarefaction waves reach the boundaries, at t — 0.2 cos a. 

Fig |A.4l shows the primitive variable profiles for all schemes against a one-dimensional reference solution obtained 
on a base grid of 1024 zones with 5 levels of refinement. Errors in L\ norm, computed with respect to the same 
reference solution, are sorted in Table [2] for density and the normal component of magnetic field. The out-coming 
wave pattern is comprised, from left to right, of a fast rarefaction, a compound wave (an intermediate shock followed 
by a slow rarefaction), a contact discontinuity, a slow shock and a fast rarefaction wave. We see that all discontinuities 
are captured correctly and the overall behavior matches the reference solution very well. The normal component of 
magnetic field is best described with MP5 and does not show erroneous jumps. Indeed, the profiles are essentially 
constant with small amplitude oscillations showing a relative peak ~ 0.7%. Divergence errors, typically < 10~ 2 , 
remain bounded with resolution and tend to saturate when the damping parameter On > 0.4 for both 2 and 3D 
calculations, see Fig |A.7l In this sense, our results favourably compare to those of HQ and & 

Fifth-order methods exhibit less dissipation across jumps, with fewer points in each discontinuous layer. Still, the 
accuracy gained from third to fifth-order accurate schemes (see Table |2]l is only a factor 1 .5 - 2 since interpolation 
across discontinuities usually degenerates to lower-order to suppress spurious oscillations. 



4.2.2. Three-Dimensional Shock tube 

The second Riemann problem was introduced by [43] and later considered by 1 44 , 5210] and by J28, 39] in 3D. 
The primitive variables are initialized as 



V, = 



1.08,1.2,0.01,0.5, 
2 4 



3.6 



1,0,0,0, 



V47T V47T V47T 



V47T V47T V47T 

T 

I 



,0.95 for X!<0 



(50) 



for xi > 



where V = (p, Vj, V2, v^, Bi, B2, B3, p). A reference solution at t — 0.2 is obtained on the domain x € [-0.75,0.75] 
using 2048 grid points and 5 levels of refinement. Our setup draws on the three dimensional version of 112 811 and 
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M39I1 where the initial condition (T50b is rotated using Eq. d45b by the angles a and y such that tancr — -1/2 and 
tan y = 1/(2 v5) (corresponding to tan/? = 1 /4). With this choice the planar symmetry is respected by an integer shift 
of cells. The computational domain consists of 768 x 8 x 8 zones and spans [-0.75,0.75] in the x direction while 
y,z & [0, 0.015625]. Computations stop at t — 0.2 cos a cosy (note the misprint in ll39ll ). 

Fig |A.5l and lA.6l show primitive variable profiles obtained with third- and fifth-order schemes, respectively. The 
wave pattern consists of a contact discontinuity that separates two fast shocks, two slow shocks and a pair of rota- 
tional discontinuities. Table [2] confirms again that the gain from high-order methods is not particularly significant 
when the flow is discontinuous. Our results favorably compare with those of other investigators and no prominent 
over/under-shoots are observed. Moreover, the amount of oscillations in the normal component of the magnetic field 



is comparable to (or smaller than) those found in 128L 13911 and divergence errors behave in a very similar way to the 
2D case (see also the right panel in Fig |A.7l ). 

The computational costs relative to that of WENO+3 (= 1) are found, for this problem, to be 0.99 : 1.48 : 1.25 
forLim03 , WENO-Z and MP5, respectively. 

4.3. Iso-density MHD Vortex advection 

The following problem has been introduced in f 5j] and lately considered by 0. l2lll . The initial condition, 
satisfying the time-independent MHD equations, consists of a magnetized vortex structure in force equilibrium that 
propagates along the main diagonal of the computational box (a square in 2D and a cube in 3D). Here we set a p = 0.4. 

4.3.1. Two Dimensional Propagation 



Following Dumbser et al. 12111 . we perform computations on the Cartesian box [-5,5] 2 with an initial flow 



describedbyp = l,v = l+(-y, x, 0) Ke q{l ~ r \ B = (-y, x, 0)//e 9<1 -' _) and p = 1 + 1 /(4 q) (p 2 (1-2 q r 2 )-K 2 p) e 2 « (1 -''~>. 
The constants k and p are chosen to be equal to 1/27T while r = sjx 2 +y 2 . The simulations are evolved for 10 time 
units with periodic boundary conditions, i.e. a single passage of the vortex through the domain. The parameter q is 
chosen equal to 0.5 for third-order schemes, effectively reproducing the configuration shown in JHQ]. For WENO-Z 
and MP5, on the other hand, we choose q = 1 in order to reduce the unwanted effects produced by the small jump in 



the magnetic field at the periodic boundaries, as argued in |2 111 . 



In order to compare our results to the findings of the latter study, we report, in Table[3] errors for B x measured both 
in L\ and Ln norms and the corresponding convergence rates. All schemes quickly converge to the asymptotic order 
of accuracy. Remarkably, errors obtained with the third-order schemes are identical and somewhat better than those 
of 0]. At the resolution of 128 2 , fifth-order schemes yield errors ~ 4 times smaller than third-order ones at 256 2 . A 
comparison between third- and fifth-order schemes from Fig |A.12l reveals that divergence errors rapidly decrease with 
resolution following a similar pattern. This eloquently advocates towards the use of higher-order schemes. 

4.3.2. Three Dimensional Propagation 

We propose a novel three dimensional extension of the vortex problem, consisting of similar initial conditions as 
the 2D case, albeit the radius r now refers to the spherical one, r = yjx 2 + y 2 + z 2 . The perturbation of pressure is now 
given by 



4q 



p 2 (l-2q(r 2 -z 2 ))-K 2 p 



e 2?(1 - r2) , (51) 



while we prescribe also a vertical velocity v z = 2. The computational domain is the cube [-5,5] 3 with periodic 
boundary conditions. The evolution stops after 10 time units. 

The last four columns of Table [3] report the L\ and L2 norm errors of B x showing an excellent agreement with 
the analytical solution. Notice that the errors measured in L2 norm are systematically smaller than L\ errors and a 



comparison between similar configurations using different norms (as reported in 112 111 ) may be deceitful. Keeping that 
in mind and given the somewhat diverse configurations, one can see that our results (in Ln norm) are competitive with 
those of 12111 at least at a qualitative level. 

Divergence errors, shown in Fig |A.12| quickly decrease as the mesh thickens and fall below 10~ 8 at the resolution 
of 128 3 for the fifth-order schemes. 

The computational cost is in accordance with previous tests, giving a ratio of 1 : 0.99 : 1 .46 : 1 .24 for WENO+3, 
Lim03 , WENO-Z and MP5, respectively. 
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Table 3: L\ and L2 norm errors and corresponding convergence rates for the MHD Vortex problem in 2D (columns 3-6) and 3D (columns 7-10) at 

t = 10. 

Two Dimensions Three Dimensions 



Method 


N x 


ei(B x ) 






Ol 2 


ei(B x ) 


U 


e 2 (B x ) 


Ol 2 


WENO+3 


32 


2.49E-03 


- 


1.94E-04 


- 


7.81E-04 


- 


1.74E-05 


- 




64 


4.13E-04 


2.6 


1.73E-05 


3.5 


1.29E-04 


2.6 


1.12E-06 


4.0 




128 


5.72E-05 


2.9 


1.16E-06 


3.9 


1.82E-05 


2.8 


5.38E-08 


4.4 




256 


7.69E-06 


2.9 


7.24E-08 


4.0 


- 


- 


- 


- 


Lim03 


32 


2.49E-03 


- 


1.94E-04 


- 


7.81E-04 


- 


1.74E-05 


- 




64 


4.13E-04 


2.6 


1.73E-05 


3.5 


1.29E-04 


2.6 


1.12E-06 


4.0 




128 


5.72E-05 


2.9 


1.16E-06 


3.9 


1.82E-05 


2.8 


5.38E-08 


4.4 




256 


7.69E-06 


2.9 


7.24E-08 


4.0 










WENO-Z 


32 


8.17E-04 




1.02E-04 




1.63E-04 




7.39E-06 






64 


5.10E-05 


4.0 


2.89E-06 


5.1 


1.07E-05 


3.9 


1.50E-07 


5.6 




128 


1.83E-06 


4.8 


5.23E-08 


5.8 


3.78E-07 


4.8 


1.87E-09 


6.3 




256 


5.94E-08 


4.9 


8.28E-10 


6.0 










MP5 


32 


9.57E-04 




1.04E-04 




1.96E-04 




7.34E-06 






64 


5.16E-05 


4.2 


3.02E-06 


5.1 


1.07E-05 


4.2 


1.53E-07 


5.6 




128 


1.75E-06 


4.9 


5.15E-08 


5.9 


3.66E-07 


4.9 


1.85E-09 


6.4 




256 


5.69E-08 


4.9 


8.04E-10 


6.0 











4.4. Advection of a magnetic field loop 

We now consider the advection of a magnetic field loop. For sufficiently large plasma B, specifying a thermal 
pressure dominance, the loop is transported as a passive scalar. The preservation of the initial circular shape tests the 
scheme's dissipative properties and the correct discretization balance of multidimensional terms [27, 28, 32, 39|. 



4.4.1. Two Dimensional Propagation 

Following Il27l,l22ll . the computational box is defined by x e [-1,1] and y e [-0.5, 0.5] discretized on 2N V x N y 
grid cells (N y = 64). Density and pressure are initially constant and equal to 1. The velocity of the flow is given 
by v = Vo(cosa, sin a) with Vq = V5, sin a = 1/ Viand cos a- = 2/ y5. The magnetic field is defined through its 
magnetic vector potential as 



A 7 = 



ao + a 2 r 
A (R - r) 




if < r < Ri , 
if Ri < r < R, 
if r > R , 



(52) 



where A = 10 -3 , R = 0.3, Ri = 0.2R, a 2 = -0.5A /R U a = A (R - Ri) - a 2 R\ and r 



yjx 2 + y 2 . The modification 



to the vector potential in the r < Ri region (with respect to similar setups presented by other investigators) is done 
to remove the singularity in the loop's center that can cause spurious oscillations and erroneous evaluations of the 
magnetic energy. The simulations are allowed to evolve until t — 2 ensuring the crossing of the loop twice through 
the periodic boundaries. 

In Fig. IA.8l the magnetic energy density is displayed for the Lim03 , WENO+3, WENO-Z and MP5 schemes, 
along with iso-contours of the z component of the magnetic vector potential. The initial circular shape is preserved 
well by all schemes. The third-order schemes are substantially more diffusive, as can be seen on the borders of the 
loop. This is confirmed by the time evolution of the magnetic energy density (normalized to its initial value), plotted 
in the left panel of Fig lA.9l The power law behaviour is similar for the schemes of the same order, with the MP5 
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method being the least diffusive. No pronounced difference is found between the Lim03 and WENO+3 schemes, for 
this particular problem. 

The divergence of magnetic field measured in L\ norm is shown in the left panel of Fig IA.11I as a function of 
a p e [0, 1 .0]. For the fifth-order schemes errors are minimized when a p > 0.4 whereas Lim03 and WENO+3 present 
smaller errors for a p < 0.2. 

4.4.2. Three Dimensional Propagation 

The three dimensional version of this problem is particularly challenging as the correct evolution depends on how 
accurately the V ■ B = condition is preserved and how the multidimensional MHD terms are balanced out. The 
computational domain -0.5 < x < 0.5, -0.5 < y < 0.5, -1.0 < z < 1.0 is resolved onto 128 x 128 x 256 zones. As for 
the two-dimensional case the vector potential A3 is used to initialize the magnetic field, which is then rotated using 
the coordinate transformation given by Eq. d45b with a = and y = tan -1 1 /2. Even though the loop is rotated only 
around one axis, the velocity profile (v A , v y , v z ) = (1,1,2) makes the test intrinsically three-dimensional. Once again, 
pressure and density are taken uniform and equal to unity while boundary conditions are periodic in all directions. 

The preservation of the loop's shape can be seen in Fig. IA. 101 All schemes preserve the shape, with Lim03 and 
WENO+3 being equally more diffusive (notice the thickness of the dark area at the loop's borders, as well as the 
brighter ring just inside the loop). As for the 2D case, one can see that MP5 is the least diffusive in preserving the 
magnetic energy (right panel of Fig. IA.9I ), while the dissipation rates for Lim03 and WENO+3 practically coincide. 
Moreover, the three-dimensional L\ norm error of V • B (right panel of Fig IA.1U exhibits a behaviour similar to the 
two dimensional case. As before, the relative CPU scaling between WENO+3, Lim03 , WENO-Z and MP5 for this 
test problem is 1 : 0.97 : 1.45 : 1.25. 



4.5. Orszag-Tang 

The Orszag-Tang vortex system describes a doubly periodic fluid configuration leading to two-dimensional su- 
personic MHD turbulence. The domain [0, l] 2 is initially filled with constant density and pressure respectively 
equal to p — T 2 and p = F, while velocity and magnetic field are initialized to v = (- sin 2jry, sin 2nx, 0) and 
B = (- sin27ry, sin47o:, 0), respectively. Although an analytical solution is not known, its simple and reproducible set 
of initial conditions has made it a widespread benchmark for inter-scheme comparison, see for example ll52ll . Density 
contour plots, as in ll32ll are shown in the top and bottom rows of Fig. IA.13l at t = 0.5 and t = 1, respectively, using a 
resolution of 256 2 points. The dynamics is regulated by multiple shock interactions leading to the formation of small 
scale vortices and density fluctuations. Our results at t = 0.5 are in good agreement with previous investigations, e.g. 
1 3oL 52, 35, 42, 32], with WENO+3 and Lim03 showing increased numerical dissipation when compared to WENO- 
Z and MP5. This is further confirmed in Fig |A.15l where horizontal cuts at y = 0.3125 in the pressure distribution are 
plotted against a reference solution obtained with the second-order CT-CTU scheme of j27[] on a finer mesh (1024 2 ), 
see also OMEHS- 



The most noticeable difference occurs at t = 1, when the fifth-order schemes (in particular, MP5) reveal the 
formation of a central magnetic island featuring a high density spot also recognizable in the results of 13211 and 
in l38ll for the isothermal case. This structure is absent in the third-order schemes and may be induced by the 
decreased effective resistivity across the central current sheet, as discussed in lT3Sh - 

Divergence errors, shown in Fig IA. 141 at t — 0.5, are comparable with those given by other investigators (e.g. 
|42j,|33[]) and reach their maximum magnitude in presence of discontinuous features. 

The computational cost of Lim03 , WENO-Z and MP5 relative to that of WENO+3 (= 1) are found to be 
1.01 : 1.47 : 1.29, in analogy with the previous results. 



4.6. Kelvin-Helmholtz Unstable Flows 



As a final example, we propose the nonlinear evolution of the Kelvin-Helmholtz instability in two dimensions. 
The base flow consists of a single shear layer with an initially uniform magnetic field lying in the xz plane at an angle 
6 = 7r/3 with the direction of propagation: 



^tanh(4o,0 

2 \yo 



B = C a ^P 



cos 6, 0, sin 



(53) 
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where M — 1 is the Mach number, yo — 1/20 is the steepness of the shear, c a — 0.1 is the Alfven speed. Density and 
pressure are initially constant and equal to p — 1 and p — 1 /T. A single-mode perturbation v y = v v o sin (2nx) exp [-y 2 /cr 2 J 
with Vyo = 10 2 , cr = 0.1 is super-imposed as in 0611 . Computations are carried out in a Cartesian box [0, 1] x [-1, 1] 
for t = 20 time units on a N x X 2N X mesh, where N x = 64, 128, 256. 

The evolutionary stages are shown in Fig IA.161 where we display color maps of the ratio {B 2 + B 2 )'- /B- at the 
largest resolution 256 x 512 for WENO+3, Lim03 , WENO-Z and MP5. For t < 5 the perturbation follows a linear 
grow th phase during which magnetic field lines wound up through the formation of a typical cat's eye vortex structure, 
1 36l 31], see the top row in Fig lA. 161 During this phase, magnetic field lines become distorted all the way down to the 
smaller diffusive scales and the resulting field amplification becomes larger for higher magnetic Reynolds numbers. 
As such, we observe in the top row of Fig IA. 171 that the magnetic energy grows faster not only as the resolution is 
increased from 64 to 256 mesh points (green, red, black), but also when switching from a third-order to a fifth order 
scheme (solid vs. dotted lines). In particular, one can see that half of the grid resolution is needed by MP5 to match 
the results obtained with WENO+3. A somewhat lesser gain can be inferred by comparing WENO-Z and Lim03 . 
Similarly, the growth rate (computed as Av v = (v^ - v 5 " ■ )/2 see bottom panel in Fig. IA. 17b . is closely related to 
the poloidal field amplification and evolves faster for smaller numerical resistivity and thus for finer grids and/or less 
dissipative schemes. 

Field amplification is eventually prevented when t > 8 by tearing mode instabilities leading to reconnection events 
capable of expelling magnetic flux from the vortex (second row in Fig. IA.161 I. [31]. Throughout the saturation phase 
(third and fourth row in Fig lA. 16b the mixing layer enlarges and the field lines thicken into filamentary structures. 
During this phase one can clearly recognize that small scale structures are best spotted with the fifth-order methods 
while they appear to be more diffused with WENO+3 and Lim03 . 

The CPU costs relative to that of WENO+3 (= 1) follow the ratios 0.98 : 1 .48 : 1 .24 for Lim03 , WENO-Z and 
MP5, respectively, and confirm the same trend already established in previous tests. 



5. Conclusions 

We have presented a class of high-order finite difference schemes for the solution of the compressible ideal MHD 
equations in multiple spatial dimensions. The numerical framework adopts a point-wise, cell centered representation 
of the primary flow variables and has been conveniently cast in conservation form by providing highly accurate 
interface values through a one-dimensional finite volume reconstruction approach. The divergence-free condition of 
magnetic field is monitored by introducing a scalar generalized Lagrange multiplier, as in M20I1 . offering propagation 
as well as damping of divergence errors in a mixed hyperbolic/parabolic way. This greatly simplifies the task of 
obtaining highly accurate solutions since the reconstruction process can be carried out on one-dimensional stencils 
using the information available at cell centers. In this respect, our formulation completely avoids expensive elliptic 
cleaning steps, does not require genuinely multidimensional interpolation and eludes the complexities required by 
staggered mesh algorithms. Selected numerical schemes based on third- as well as fifth-order accurate constraints 
have been presented and compared. 



The recently improved version of the third-order WENO scheme (WENO+3, 15311 ') and the Lim03 reconstruction 
based on new limiter functions (introduced in OH) perform equally well exhibiting third-order accuracy in 
smooth problems and non-oscillatory transitions at discontinuities. 

The new fifth-order WENO scheme (WENO-Z, see |0|) and the monotonicity preserving algorithm (MP5) of 



1 48|] yield high-quality results on all of the selected tests and report orders of accuracy close to 5 for multidi- 
mensional smooth problems. Both WENO-Z and MP5 perform with a greatly reduced amount of numerical 
dissipation and provide highly accurate solution with much fewer grid points when compared to third-order 
accurate schemes. Still, we have found MP5 to give slightly better results WENO-Z in terms of reduced com- 
putational cost, improved accuracy and sharper transitions at discontinuous fronts. 

Fifth-order schemes are found to be < 50 (for WENO-Z) and < 30 (for MP5) per cent slower than third-order 
ones, depending on the particular choice. This favorably advocates towards the use of higher order schemes 
rather than lower order ones, since the same level of accuracy can be attained at a much lower resolution still 
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giving a tremendous gain in computing time. For three-dimensional problems, for example, the gain can be 
almost two orders of magnitude in CPU cost. 



The results obtained with the present finite difference formulation are competitive (in terms of accuracy and 
description of discontinuities) with recently developed FV schemes (e.g., J2II and noticeably improve over 
traditional 2 nd order Godunov-type schemes in terms of reduced numerical dissipation. The benefits offered by a 
high-order method such as the ones presented here are particularly relevant in the context of MHD applications 
involving both smooth and discontinuous flows. 
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Appendix A. Conservative eigenvectors of the GLM-MHD Equations 



The 9x9 matrix of the conservative MHD equations in one dimension introduced can be decomposed, given the 
eigenvalues (see Eq. 13, to the corresponding left and right eigenvectors. Following partially the notation of 1 41 , 30ll . 
we define 



2 2 
a - ci 



a f = 



2 2 

2 c f~ a ~ 
a, = 



2 ? 
c — c 



By = 



B 



(A.l) 



where a = ^Tplp denotes the speed of sound. With this notation, the right eigenvectors in matrix form will be given 
by 
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where S = sign(B x ), H ftS = a f , s (0.5v 2 + c 2 f s - y 2 a 2 ), J ftS0 = a s jc sJ S and J ftS \ = a s jap? . 
On the other hand, the left eigenvectors are given by 
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1-1,9 = 



2c h 
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where we prescribe t = (V - I) I a 2 , j\ = (T- l)/2, 72 = (T-2)/(T- 1) and /(/,.,) („ f ,B f ) = F 1 a^ s (vi,B i ), with i = 
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One Dimension 



Three Dimensions 
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Figure A. 1: L\ norm errors computed for the one-dimensional Alfven wave propagation (left panel) and the rotated three-dimensional version 
(right panel). The cross, triangle, plus sign and square symbols refer to computations carried out with WENO+3, Lim03 , WENO-Z and MP5, 
respectively, at the resolution 16, 32, 64, 128 and 256 points using a CFL number of 0.8 (in ID) and 0.3 (in 3D). The dotted lines gives the ideal 
convergence slope, that is, oc A.v 3 and oc A* 5 , respectively. 
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Figure A. 2: Scatter plots of the y component of magnetic field in the original one-dimensional frame at r = 5/3, after 5 revolutions. Each panel 
plots every point of the three-dimensional array — B x sin a + B, cos a as a function of the longitudinal coordinate k ■ x/|k| along the direction of wave 
propagation. The lack of scatter demonstrates that the algorithm retains the expected planar symmetry. The solid line gives the reference solution 
at t = while dotted, dashed and dot-dashed lines corresponds to computations carried with N x = 16, 32, 64 points, respectively. The CFL number 
was set to C a = 0.3. 
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Figure A. 3: Long term decay of circularly polarized Alfven waves after 16.5 time units, corresponding to ~ 100 wave periods. In the left panel, 
we plot the maximum value of the vertical component of velocity as a function of time for the WENO-Z (solid line) and WENO+3 (dashed line) 
schemes. For comparison, the dotted line gives the result obtained by a second-order TVD scheme. The panel on the right shows the analogous 
behavior of the vertical component of magnetic field B. for Lim03 and MP5. For all cases, the resolution was set to 120 X 20 and the Courant 
number is 0.4. 
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Figure A. 4: Primitive variable profiles for the 2D shock tube problem at t = 0.2 cos a = 0.2/ V2, along the rotated direction xi . From left to right: 
density, transverse velocity, longitudinal and transverse magnetic field components are displayed. The mesh resolution is 600 X 6 and the Courant 
number is 0.4. Symbols correspond to the 2D computations whereas the solid lines gives the reference solution. 
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Figure A. 5: 


Primitive variable profiles for the 3D shock tube problem at t = 
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Density, pressure, velocity and magnetic field components parallel and transverse to the direction of propagation are plotted as functions of the 
longitudinal component x\ . The mesh resolution is 768 X 8 X 8 and the Courant number is 0.3. 
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Figure A. 6: Same as Fig |A.5l but for the fifth-order schemes WENO-Z and MP5. 
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Figure A. 7: Divergence errors as function of the damping parameter a p for the shock tube problems in 2D (left, t = 0.2/ V2) and 3D (right, 
t = 0.8/ V2T). Symbols in black color are used to distinguish between different schemes at the nominal resolutions (600 X 6 in 2D and 768 X 8 X 8 
in 3D), see the legend. Computations carried at twice the resolution (1200 X 12 in 2D and 1536 X 16 X 16 in 3D) are shown using symbols in red 
color. 




Figure A. 8: Magnetic energy density for the 2D field loop problem at t = 2 computed with the third-order (left) and fifth-order (right) schemes at 
the resolution of 128 X 64 points with Courant number C a = 0.4. Magnetic field lines are overplotted using 9 contour levels equally spaced between 
10~ 5 and 10~ 3 . 
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Figure A. 9: Time evolution of the magnetic energy density, normalized to its initial value, for the 2D (left) and 3D (right) field loop problem at 
the resolution of 128 X 64 grid points. The magnetic energy is better conserved for the MP5 method. Lim03 and WENO+3 show no pronounced 
difference for this particular problem. 




Figure A. 10: Magnetic energy density for the 3D field loop problem at t = 1 computed on 128 X 128 X 256 grid zones with Courant number 0.3. 
From left to right: Lim03 , WENO+3, WENO-Z, MP5. All schemes preserve the circularity of the loop, with the fifth-order schemes displaying 
sharper borders. 
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Figure A. 1 1 : Divergence errors as function of the damping parameter a p for the field loop test problem in 2D (left, t = 2) on 1 28 X 64 grid points 
and 3D (right, t = 1) on 128 X 128 X 256 grid points. 
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Figure A. 12: L\ norm error of the divergence of magnetic field as functions of the resolution (N x ) for the 2D (left panel) and 3D (right panel) vortex 
problems at t = 10. Different symbols corresponds to the selected reconstruction algorithms. 
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Figure A. 13: Density contour plots for the Orszag-Tang system at ; = 0.5 (top) and t = I (bottom) for the selected schemes using 256 2 grid points. 
Thirty equally spaced levels ranging from 0.383ir 2 to 2.2414r 2 for the top panel and from 0.1944r 2 to 1.9337r 2 for the bottom panel are shown. 
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Figure A. 14: Divergence errors for the four selected scheme att = 0.5 on 256 2 grid zones. 
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Figure A. 16: Snapshots of the evolution of the Kelvin-Hemlholtz unstable layer at t = 5 (first panel from top), t = 8 (second panel), t = 12 (third 
panel) and t = 20 (bottom panel). The images show the ratio of the poloidal field strength and the toroidal component, -^B]. + By/B z . Left to right 
columns corresponds to computations obtained with WENO+3, Lim03 , WENO-Z and MP5, respectively, at the resolution of 256 X 512. Note 
how the colorbar maximum value changes at different instant to reflect tj^ corresponding magnetic field strength. 
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Figure A. 17: Volume integrated magnetic energy (top panels) and growth rate (computed as Av v = (v y m[ 



)/2) as functions of time. Here 



BZ = B~ + B- accounts for the "poloidal" contribution only. Solid and dotted lines corresponds to integrations carried with WENO-Z and 
Lim03 (left panels), MP5 and WENO+3 (right panels). The different colors, green, red and black indicate different numerical resolution, i.e., 64, 
128 and 256, respectively. 
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